Particle Filter with SynthTown data

In this example, we track the numbers of vehicles on different links in a transportation network from the observations about probe vehicles. We have synthesized the following dataset from MATSIM equil example:

The transportation network looks like the following. Individuals in this scenario leave their homes on link 1 to their workspace on link 20 around 6:30am, and leave their workspace to their home around 7pm. Our goal is to track the numbers of vehicles at different locations in a transportation network from the numbers of probe vehicles at link 1 and link 20 using a particle filter.

synthesized data from MATSIM equil

The particle filter algorithm takes 3 pieces of information as input: - h(x,t) the event rates of the uniformized stochastic kinetic process, - S.uniformized the stoichiometric matrix, and - obs.prob(x,t) the likelihood of latent state with respect to observations.

It generates as output: - Xt_est the estimated numbers of vehicles at different locations.

while(t < nrow(Xt_est)) {
  gamma = max(apply(Xt, 1, function(x) {sum(h(x,t))*gamma.inc} ))
  event.time.sched <- t+rexp(1,rate = gamma)

  # resampling/updating
  if ((trunc(event.time.sched)-trunc(t)) > 0) {
    for (i in c(1:(trunc(event.time.sched)-trunc(t)))) {
      W = obs.prob(Xt,t,i)
      Xt = Xt[sample(length(W),particles,replace=TRUE,prob=W),]
      Xt_est[trunc(t)+i,] = colMeans(Xt)
    }
  }
  
  # sampling/prediction
  Xt = t(apply(Xt, 1, function(x) {
    event.rate = c(h(x,t)/gamma, 1-sum(h(x,t))/gamma)
    vt <- sample(c(1:length(event.rate)),1,prob = event.rate)  
    x <- x + S.uniformized[,vt]
  }))
  
  # updating time
  t = event.time.sched
}

# plot result
layout(matrix(1:ncol(Xt_real),ncol=1), heights=apply(Xt_real,2,max))
par(mar=c(0,0,0,0),oma=c(5,4,4,2)+.1)
for(i in 1:ncol(Xt_real)){
  plot(Xt_real[,i],type='l',col=1,xaxt='n',xlab='',ylab=colnames(Xt_real)[i],lty=2)
  lines(Xt_est[,i], col=2,lty=1)
  if(i==1) legend('bottomright',col=1:2,lty=2:1,legend=c('groundtruth','particle filter') )
}
axis(side=1)